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Subspace Leakage Analysis and Improved 
DOA Estimation with Small Sample Size 
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Abstract 

Classical methods of DOA estimation such as the MUSIC algorithm are based on estimating the 
signal and noise subspaces from the sample covariance matrix. For a small number of samples, such 
methods are exposed to performance breakdown, as the sample covariance matrix can largely deviate 
from the true covariance matrix. In this paper, the problem of DOA estimation performance breakdown 
is investigated. We consider the structure of the sample covariance matrix and the dynamics of the root- 
MUSIC algorithm. The performance breakdown in the threshold region is associated with the subspace 
leakage where some portion of the true signal subspace resides in the estimated noise subspace. In this 
paper, the subspace leakage is theoretically derived. We also propose a two-step method which improves 
the performance by modifying the sample covariance matrix such that the amount of the subspace leakage 
is reduced. Furthermore, we introduce a phenomenon named as root-swap which occurs in the root- 
MUSIC algorithm in the low sample size region and degrades the performance of the DOA estimation. 
A new method is then proposed to alleviate this problem. Numerical examples and simulation results 
are given for uncorrelated and correlated sources to illustrate the improvement achieved by the proposed 
methods. Moreover, the proposed algorithms are combined with the pseudo-noise resampling method 
to further improve the performance. 
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I. Introduction 

Classical parameter estimation methods of direetion-of-arrival (DOA), Doppler shifts, frequen- 
eies, time delays, etc. such as the multiple signal classifieation (MUSIC) |[I1, root-MUSIC (III, 
and estimation of signal parameters via rotational invarianee teehniques (ESPRIT) [[3l are based 
on estimating the signal and noise subspaees from the sample data eovarianee matrix. It is well- 
known that these methods suffer from performanee breakdown for a small number of samples 
or low signal-to-noise ratio (SNR) values where the expected estimation error departs from the 
Cramer-Rao bound (CRB) [|4l. The SNR region at whieh this phenomenon happens is known as 
the threshold region. 

The fidelity of the sample data eovarianee matrix to the true data eovarianee matrix plays a 
critieal role in a sueeessful estimation. At the low SNR and/or small sample size region, the 
sample data eovarianee matrix ean largely deviate from the true one. There are various methods 
introduced in the literature whieh target at improving the estimation of the eovarianee matrix 

nai-iiia. 

Diagonal loading [[51 and shrinkage-based flU methods improve the estimate of the data 
eovarianee matrix by sealing and shifting the eigenvalues of the sample data eovarianee matrix. 
However, the eigenveetors are kept unehanged. As a result, the estimated signal and noise 
projection matrices from the improved eovarianee matriees are exaetly the same as those obtained 
from the sample data covariance matrix. Therefore, these methods are not really benefieial for 
the subspaee-based parameter estimation algorithms. 

Data eovarianee matrix estimation can be also improved by the means of using forward- 
backward averaging (FB) 0 and spatial smoothing-based teehniques (jH. The effect of FB is 
known to be equivalent to approximately doubling the number of samples. Thus, the eovarianee 
estimate improves accordingly. The spatial smoothing technique can also be interpreted as 
virtually increasing the number of samples at the eost of averaging over sub-arrays of smaller size 
eompared to the whole array. These approaehes ean also deeorrelate pairs (in ease of FB) or more 
eorrelated source signals. In dUl, techniques from random matrix theory have been developed 
to improve the performanee of the MUSIC algorithm. The introdueed method eonsiders the 
asymptotic situation when both the sample size and the number of array elements tend to infinity 
at the same rate. It is then inferred that the improved method gives a more accurate deseription 
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of the situation when these two quantities are finite and eomparable in magnitude. However, the 
performanee of the introdueed method is not satisfaetory at the small sample size seenario [[T3ll . 

A more promising approaeh to remedy the performanee breakdown at the threshold region was 
introdueed in IfTOl and has been further improved in [fTTlI and ITT^ . These methods are based on 
a teehnique ealled pseudo-noise resampling whieh uses synthetieally generated pseudo-noise to 
perturb the original noise. The pseudo-noise is added to the observed data, and a new estimate of 
the eovarianee matrix is obtained, whieh leads to new DOA estimates. This proeess is repeated 
for a number of times, and the final DOAs are determined based on the bank of the DOA 
estimates. 

In this paper, we taekle the problem of the performanee breakdown at the threshold region 
by eonsidering the strueture of the sample data eovarianee matrix and the dynamies of the root- 
MUSIC algorithm. It is shown in [fT4ll that the performanee breakdown problem is assoeiated 
with the inter-subspaee leakage “whereby a small portion of the true signal eigenveetor resides 
in the sample noise subspaee (and viee versa)”. The notion of leakage eomes originally from 
the performanee assessment strategy based on the first order approximation of the estimation 
error eaused by the perturbed subspaee estimate, whieh happens beeause of the additive noise 
eontribution ifTSl - lfT^ . This approaeh direetly models the leakage of the noise subspaee into 
signal subspaee and allows to eompute the eorresponding perturbation matrix between the 
eomponents of the subspaees. Here, we formally define the subspace leakage notion as a 
Frobenius norm of the perturbation matrix, and we present its theoretieal derivation. We propose 
a two-step method whieh improves the performanee of the root-MUSIC algorithm by modifying 
the sample data eovarianee matrix sueh that the amount of the subspaee leakage is redueed. 
Furthermore, we introduee a phenomenon named as root-swap whieh oeeurs in the root-MUSIC 
algorithm at the threshold region and degrades the performanee of the parameter estimation. A 
new method is then proposed to alleviate this problem. 

It will be shown that there are undesirable by-produets in the sample data eovarianee matrix 
that tend to zero as the number of samples goes to infinity. However, for a limited number of 
samples, these terms ean have signifieant values leading to a large amount of subspaee leakage. 
One possible approaeh to remedy the effeet of the undesirable eomponents is to eonsider the 
eigenvalue perturbation eaused by these terms. The ineorporation of this knowledge into the 
estimation method ean result in better estimates of the signal and noise subspaees. In this paper. 
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we propose a two-step algorithm in order to reduee the effeet of the undesirable terms. The 
introdueed method is based on estimating the parameters at the first step and modifying the 
eovarianee matrix using the estimated parameters at the seeond step. We will theoretieally derive 
the subspaee leakage at both steps. Then, it will be shown using numerieal examples that the 
subspaee leakage is redueed at the seeond step leading to better performanee. 

In the root-MUSIC method, the estimation error of the roots has a varianee whieh is pro¬ 
portional to the varianee of noise over the number of samples lfT 9 l . Therefore, at the threshold 
region, the varianee of the estimation error ean have a signifieant value whieh in turn ean result in 
a swap between a root eorresponding to a signal souree with another root whieh is not assoeiated 
with any signal souree. We dub this phenomenon as root-swap. Then, a new method is proposed 
to remedy this problem. The introdueed method eonsiders different eombinations of the roots as 
the eandidates for the signal sourees. These eandidates are then evaluated using the stoehastie 
maximum likelihood (SML) funetion, and the eombination that minimizes the objeetive funetion 
is pieked up for the parameter estimates. 

The rest of the paper is organized as follows. The system model is given and the root- 
MUSIC algorithm is briefly reviewed in Seetion UIl The two-step and root-swap algorithms are 
proposed in Seetion Hill Subspaee leakage is defined and theoretieally derived in Seetion |Wl 
Numerieal examples and simulation results are given in Seetion |Vl Seetion |Vl] eoneludes the 
paper. Appendix |A] gives an approximation for the probability of root-swap, and finally, the 
details of the subspaee leakage derivation for the two-step root-MUSIC algorithm are presented 
in Appendiees |B] and O 


II. System Model and Background 


An example of a noise-eorrupted linear superposition of K undamped exponentials reeeived by 
M (M > K) antennas is the array proeessing model. Thus, eonsider K number of narrowband 
plane waves impinging on a uniform linear array (ULA) from direetions 6*i, 62, - ■ ■ , Ok- Without 
loss of generality, assume —7r/2 < 6*1 < 6*2 < < Ok < t^/‘ 2 -. The antenna elements are 


separated from eaeh other by a distanee of d < A/2 where A is the wavelength of the plane 


wave impinging on the array. The steering veetor of the array a{ 0 ) G given as 



( 1 ) 
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where (■)^ stands for the transposition operator. At time instant t E N, the reeeived veetor 
x{t) e is given by 

K 

x{t) = ^ a{ei)si{t) + n{t) (2) 

i=l 

where Si{t) G C is the amplitude of the i-\h wave (souree) and n{t) G is the noise veetor at 
time t. By arranging the amplitudes of the sourees in the veetor s{t) = [si(f), S2(f), • • • , G 

forming the Vandermonde matrix A = [a(6'i), a(6*2),--- , a{6K)] G the 

model Q ean be rewritten in matrix-vector form as 

x{t) = As{t) + n{t). ( 3 ) 

We consider the noise vector n{t) to be independent from the sources and noise vectors at 
other time instances and to have the circularly-symmetric complex jointly-Gaussian distribution 
A/c( 0 , u^Im) where Im is the identity matrix of size M. Considering the system model dH), the 
data covariance matrix R G is given by 

E {x{t)x^{t)} = ASA^+ allM ( 4 ) 

where S = E G is the source covariance matrix and (■)^ and E{-} stand 

for the Hermitian transposition and the expectation operators, respectively. 

Let Ai < A2 < < Am be the eigenvalues of R arranged in nondecreasing order, and 

let Qi, g2, - , 9m-k the noise eigenvectors associated with Ai, A2,-- - , Xm-k and 

Cl , 62, • • ■ , bk be the signal eigenvectors corresponding to Xm-k+i, Xm-k+2, ■ " ^ Xm- 
Let also G G and E G be defined as G = [g^, 92, - ' ; 9m-k\ 

E = [ei, 62, • • • , bk]- The range spaces of G and E represent the true noise and signal 
subspaces, respectively. 

Let N number of snapshots (samples) be available. The basic method for estimating the data 
covariance matrix from the samples x{t) (1 < f < N) is 

1 ^ 

R = — ^x{t)x^{t) ( 5 ) 

' t=i 

where R G is the sample data covariance matrix. Consider the eigendecomposition of 

R. Let g^, 92 ^ ' 1 Om-k be the estimated noise eigenvectors and ei, 62, • • • , be the 
estimated signal eigenvectors. Form G G and E G by placing the estimated 
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noise and signal eigenvectors as the columns of G and E, respectively. The range spaces of G 
and E represent the estimations of the noise and signal subspaces, respectively. 

Recalling dU) and defining 2 ; = gi27r(d/A)sin(e)^ steering vector can be rewritten as a{z) = 
[ 1 , • • • , . In the root-MUSIC method, the roots of the equation aF{z~^)GG a{z) = 

0 which are located inside the unit circle are considered. These roots are sorted based on their 
distance to the unit circle, and the first K number of the roots which are closer to the unit 
circle are picked. The estimates of the DOAs denoted by 6^1, 6*2, • • • , Ok are then obtained by 
multiplying the angles of the selected roots by X/{27rd) and taking the inverse sinusoid function 
of the results. 


III. Proposed Methods 
A. Two-step root-MUSIC algorithm 

Let us start by expanding Q using ([3]) as follows 

^ 1 ^ 

i^= — ^ {As{t) + n{t)) {As{t) + n{t)f 

' t=i 

f 1 ^ 1 1 ^ 

t=i J ' t=i 

Comparing ® with dH), it can be observed that the expansion of R consists of four terms 
while the model for R comprises two summands. The first two terms of R given by db]) can 
be considered as estimates for the two summands of R, which represent the signal and noise 
components, respectively. The last two terms of i? in db]) are undesirable by-products which can 
be viewed as estimates for the correlation between the signal and noise vectors. In the system 
model under study, we consider the noise vectors to be zero-mean and also independent of the 
signal vectors. Therefore, the signal and noise components are uncorrelated to each other. As a 
result, for a large enough number of samples N, the last two terms in db]) tend to zero. However, 
the number of available samples can be limited in practical applications. In this case, the last 
two terms in (jb]) may have significant values, which causes the estimates of the signal and noise 
subspaces to deviate from the true signal and noise subspaces. 
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The main idea of our two-step root-MUSIC algorithm is to modify the sample data covariance 
matrix at the second step based on the DOA estimates obtained at the first step. The modified 
covariance matrix is obtained by deducting a scaled version of the estimated undesirable terms 
from the sample data covariance matrix. 

We derive the steps of the proposed method for a general source covariance matrix S, so 
that correlated sources can also be handled by the algorithm. Furthermore, the proposed method 
can also be beneficial in the case that the assumption on no correlation between the source and 
noise vectors is not fully met. This is achieved by estimating and removing the correlation terms 
between the source and noise vectors from the sample data covariance matrix. 

The steps of the proposed method are listed in Table HI The algorithm starts by computing the 
sample data covariance matrix ([5]). Then, DOAs are estimated using the root-MUSIC algorithm. 
The superscript refers to the estimation made at the first step. At the second step, the Van¬ 
dermonde matrix is formed using the available estimates of the DOAs. Then, the amplitudes of 
the sources are estimated such that the squared norm of the differences between the observations 
and the estimates are minimized. The corresponding problem is formulated as 


s(t) = arg min ||a;(f) — As^\ 


(7) 


The minimization of © is performed using the least squares (LS) technique and the corre¬ 
sponding solution is given as 



( 8 ) 


The noise component is then estimated as the difference between the estimated signal and the 
observation made by the array, i.e.. 


h{t) = x{t) — As(t). 


(9) 
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After estimating the signal and noise veetors, the third term in ® ean be found as 

" S (a"2) A''x(t) - x''(t)A (^a^a) a" 

_ r 1 iv ^ 

= {Im - Pa) 

t=l 

= PaRPa (10) 

where 

Pa = a(^A A) A (11) 

is an estimation for the projection matrix of the signal subspace, and 

PA ~ ~ PA ( 12 ) 

is an estimation for the projection matrix of the noise subspace. The forth term in ® is equal to 
the Hermitian of the third term, i.e., . Finally, the modified data covariance matrix is obtained 

by deducting a scaled version of the estimated terms from the initial sample data covariance 
matrix as follows 

= ^ - 7 (T + T^) . (13) 

The scaling factor 7 in (fT3l) is a real number between zero and one. Ideally, the value of 7 
would be equal to 1 if the estimates of the undesirable terms were perfect. However, estimation 
errors are inevitable, and therefore, we have introduced 7 to deal with the imperfections. The 
scaling factor 7 can be considered as a reliability factor which takes a value close to 1 for an 

estimate of T with small error and a small value if an estimate of T is erroneous. Given a value 

^( 2 ) 

for 7, the modified data covariance matrix R is computed and the DOAs are estimated again 
using the root-MUSIC algorithm. 

The value of 7 can be fixed to a predetermined value before running the algorithm, or it can 
be obtained based on the observations. Since 7 is a real number between zero and one, we can 
consider different values for 7 taken on a grid (e.g. 7 = 0, 0.1, 0.2, • • • , 1). For each value of 
7, a set of DOA estimates is obtained based on the modified data covariance matrix. Next, we 
determine which value of 7 results in a better estimation. This can be done by choosing a set 
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of DOA estimates that has a higher likelihood of being the set of true DOAs. In other words, 


we use the maximum likelihood (ML) eriterion to evaluate the quality of the estimated DOAs. 
Sinee the system model given in (jH) is stoehastie, we use the stoehastie ML (SML) objeetive 


function given by [|20ll 


FsmM = Indet 


^( 2 )-( 2 ) 

PaRPa + 


M-K 


i-L(2) 




(14) 


where Tr{-} stands for the trace operator, is an estimation of the projection matrix of the 

signal subspace obtained from the estimated DOAs based on the modified data covariance matrix 

^±( 2 ) ^( 2 ) 

and = Im — Pa ■ The objective function in (1141) is evaluated for each value of 7. Then, 

the set of DOA estimates corresponding to the value of 7 that minimizes (fT4l) is chosen as the 


output of the algorithm. 


B. Root-swap root-MUSIC algorithm 

Consider the root-MUSIC polynomial {z~^)GG^a{z) which is formed by the noise eigen¬ 
vectors obtained from the eigendecomposition of the data covariance matrix R. This polynomial 
has K number of roots on the unit circle which correspond to the signal sources. Let these K 
roots be denoted by zi, 2 : 2 , • • • > and be referred to as the true signal roots. The polynomial 
also has additional M — K — 1 number of roots inside the unit circle. Let these roots be referred 
to as the true noise roots and be denoted by zk+i, zk+ 2 , • • • , zm-i- 

An estimation for the root-MUSIC polynomial can be formed using the noise eigenvectors 
obtained from the sample data covariance matrix R. Let us assume that in the estimation of the 
noise and signal subspaces, no subspace swap has occurred BH. The estimated polynomial is 

H 

given by a^{z~^)GG a(z). This polynomial has M — 1 number of roots inside the unit circle. 
Let zi, Z 2 , - ■ ■ , 5k be the roots of the estimated root-MUSIC polynomial which correspond 
to zi, 2 : 2 , • • • , zk- We refer to these roots as the estimated signal roots. Furthermore, let 
zk+i, 5 k+2 , • ■ • ) zm -1 be the roots corresponding to zk+i, zk+ 2 ,-" : zm-i- These roots 
are referred to as the estimated noise roots. 

In the root-MUSIC method, we do not have the knowledge about which of the roots of the 
estimated root-MUSIC polynomial correspond to the true signal roots. The conventional rule is 
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TABLE I 

Two-step root-MUSIC algorithm 


Inputs: 

M, d, A, N, K, and 
received vectors ai(l), x(2), 

Outputs: 

Estimates 9 ^^, , • 


, x{N) 


9 


( 2 ) 

K 


Step 1: 

9^^\ • • • , ^ root-MUSIC ( K, d, A 

Step 2: 


a () , a [9. 


1 (1) 


, a (9 


)(i) 

'k 


A = 


Pa = a[a a) a 
Pa = — Pa 

T = PaRPa 

Determine 7 as the minimizer of (fT4l) 

,( 2 ) 


R =R--f{T + T^) 


e?, 9^^\---JP 


^root-MUSIC K ,K,d,X 


,(2) 


to select K number of the estimated roots which are closer to the unit circle as the estimates 
for the true signal roots. Then, the DOAs are estimated based on the angles of these roots. 

Due to the finiteness of the available samples, the estimated roots obtained from the sample 
data covariance matrix R deviate from their corresponding true roots obtained from the true data 
covariance matrix R. Let r, and f, represent the magnitudes of Zi and Zi for 1 < f < M — 1, 
respectively. Furthermore, let Ar* ^ fi — Vi be the difference between the magnitude of the 
f-th estimated root and the magnitude of the corresponding true root. It is shown in [fT^ that 
Arj (for the signal roots) has a variance which is proportional to cr^/N. Therefore, Ar* can 
have a significant value for a small number of samples and a large value of cx^ (low SNR 
region). Consequently, there can be a considerable probability that an estimated signal root takes 
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a smaller magnitude than an estimated noise root. We refer to this phenomenon as a root-swap. 
The root-swap probability is approximately found in Appendix as 

+ o'k\/M — K — (3/4) \ 


K M-1 


F(root-swap) ~ 1 ~ Q 


k=l m=K-\-l 




V 


(15) 


where Q {■) is the tail probability of the standard normal distribution and is the varianee 
of Arfc, and it is proportional to cr^/N. 

In the ease that the root-swap happens, seleeting the first elosest K roots to the unit eirele 
results in pieking a noise root instead of a signal root. To deal with this problem, we propose an 
algorithm that eonsiders different eombinations of the roots as eandidates for signal roots. The 
method is dubbed the root-swap root-MUSIC algorithm. 

The root-MUSIC polynomial has M — 1 number of roots inside the unit eirele. Our goal is to 
find the roots whieh have a higher likelihood of being assoeiated with the K sourees. Consider 
ehoosing K number of roots out of the M — 1 roots inside the unit eirele. There are Nc = 
(M — 1)!/ {K\{M — K — 1)!) different possible eombinations. Let T A { 0 ^^ 02 , ■ ■ ■ , 0Ar^} 
where 0i (1 < i < Nc) is a set eontaining the DOA estimates obtained from the Ath eombination 
of the roots. Then, the root-swap root-MUSIC method estimates the DOAs as 


Oi, 6 ^ 2 , • • • , 6 'x (> = arg minF sml (©) 


where Fsml (0) is the SML funetion given by 

Fsml{Q) = Indet [ PqRPq + 


Tr <! P^R 


M-K 

and P© is the signal projeetion matrix obtained from 0 as 


e 


(16) 


(17) 


P© A A(0) (A^(0)A(0)) 'a^(0). (18) 

The eomplexity of the introdueed root-swap root-MUSIC method ean be redueed by pre¬ 
eliminating some of the roots. Speeifioally, let p < K roots elosest to the unit eirele be pieked, 
and let q number of roots elosest to the origin (furthest from the unit eirele) be ignored. Our task 
is to ehoose K —p number of roots out of M—p — q — 1 roots. Then, there are = {M—p — q — 
1)!/ {{K — p)\{M — K — q — 1)!) different possible eombinations whieh is signifieantly smaller 
than Nc- The rest of the algorithm is the same as above exeept for that here eaeh eombination 
eontains K — p number of roots. Therefore, in order to evaluate the SML funetion, the fixed p 
pre-seleeted roots are added to eaeh eombination. 
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IV. Sub SPACE Leakage 

The performance breakdown of the subspace based DOA estimation methods in the threshold 
region has been associated with the subspace leakage. In this section, we study the subspace 
leakage for both steps of the proposed two-step root-MUSIC algorithm. 

A. Definition 

Recall the matrices G and E which are composed of the true noise and signal eigenvectors 
obtained from the eigendecomposition of the data covariance matrix R. Note that the matrix of 
the eigenvectors =[G E]e is a unitary matrix {QrQ^ = Im), therefore 

GG^ + EE^ = Im (19) 

or 

P^ + P = Im (20) 

where, = GG^ and P = EE^ are the true projection matrices of the noise and signal 
subspaces. 

Ideally, the estimation of each signal eigenvector ek {1 < k < K) would perfectly fall in 
the true signal subspace. In practice, however, the energy of the projection of into the noise 
subspace ||P‘*‘efc ||2 is almost surely nonzero, which can be viewed as the leakage of into the 
true noise subspace. 

We define the subspace leakage as the average value of the energy of the estimated signal 
eigenvectors leaked into the true noise subspace, i.e., 

P = (21) 

k=l 

Note that P"*" is the orthogonal projection matrix. Therefore, p can be written as 

1 ^ 

<22) 

k=l 
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Using (l20l) and some algebra, the expression (l22l) can be simplified to 

1 ^ 

P = (Im- P)h 

k=l 

■-id (!;*•■•)'■) 

= l-^Tr|pp} (23) 

^ ^- H 

where P = EE is the estimated signal projection matrix. 


B. Analysis of two-step root-MUSIC algorithm 

The estimated signal and noise projection matrices obtained from the eigendecomposition of 
the sample data covariance matrix R are deviated from the true signal and noise projection 
matrices. Let AP = P — P be the estimation error of the data covariance matrix, and let 


V ^ R- allM = ASA^ 


K 


k=l 


M-K+k 




(24) 


Denote the pseudo-inverse of V as G jj- given by 

K 




^ ^M-K+k — O'n 




(25) 


Let Pi and p 2 be the subspace leakage due to the error in the estimation of the signal and 

^ ^(2) 

noise subspaces obtained from P and P , respectively. Note that pi only depends on P and 
AP, and it is not specific to the proposed two-step root-MUSIC algorithm. 

It is shown in Appendix |B] that pi and its expected value are given by 


Pi = —TrjVfAPP^APU^} 


and 


K 


al (M - K) 
^{pi} = —^ 


A/u 


M-K+k 


NK 


2 U 


^ {^M-K+k - O-n) 


(26) 


(27) 


respectively. 
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It can be seen from (ITTI) that the expected value of the subspace leakage is proportional to 
cr^/A^. Therefore, the amount of the subspace leakage can be significant for a small number of 
samples or low SNR values. The variance of pi has also been studied in ll^ . and it has been 
shown that Var(pi) is in the order of 1/iV^. 

The subspace leakage at the second step of the two-step root-MUSIC algorithm is computed 
in Appendix O and is given by 

P 2 = (1 - 27 + 7 ^) Pi + ^ ~ Re {Tr {V^ ARP^dP} } + ^Tr {dPP^dP} (28) 

where Re{-} stands for the real part operator, and dP is the first order term in the Taylor series 
expansion of Pa around the true DOAs. It is also shown in Appendix O that the expected value 
of p 2 for a fixed value of 7 is given by 


E { 92 } = (1 - 27 + 7 ^) E {pi} 

, 2 ( 7 - 7 ^) 


NK 


R, ) W A«V^RV'ay, \ 


k=l 


2j ( 






2 2 K K 

k=i i=i I P \ [a^ P 




(29) 


where ujk — 27r((i/A) sin( 6 *fc), is a shorthand notation for a{ 6 k), and G is defined 


as 


49 A_[ 0 , 2 e-^■""^^ •••, (M- 


(30) 


It can be seen in (l29l) that for 'j = 0, E {^2} reduces to E {pi} as expected, and for 7 = 1 , 
the first two terms in (l29l) are equal to zero. 


V. Numerical Examples and Simulation Results 

In this section, the performance of the proposed two-step root-MUSIC and the root-swap 
root-MUSIC algorithms is investigated and compared with the performance of the unitary root- 
MUSIC method [l22ll and the improved unitary root-MUSIC algorithm based on pseudo-noise 
resampling iTT^ . We also consider the combination of the proposed methods with the other 
methods in order to achieve further performance improvement. Compared to the root-MUSIC 
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method, the unitary root-MUSIC algorithm has a lower eomputational eomplexity as it uses the 
eigendeeomposition of a real-valued eovarianee matrix. Furthermore, the unitary root-MUSIC al¬ 
gorithm has better performanee for the ease that the sourees are eorrelated. The improved unitary 
root-MUSIC algorithm based on pseudo-noise resampling inereases the estimator eomplexity, but 
it is advantageous in removing the outliers, whieh results in better performanee. 

We eonsider K = 2 sources impinging on an array of M = 10 antenna elements from 
directions = 35° x (vr/lSO) and 6*2 = 37° x (tt/ISO). The interelement spacing is set to 
d = \/2 and the number of snapshots is iV = 10. Each source vector s{t) is considered to be 
independent from the source vectors at other time instances and to have the circularly-symmetric 
complex jointly-Gaussian distribution A/c(0, S). The source covariance matrix S is given by 


1 r 



r 1 


where 0 < r < 1 is the correlation coefficient. The SNR is defined as SNR = lOlog^p ('^s/'^n)- 
The performance of the proposed algorithms is investigated by considering the subspace 
leakage, mean squared error (MSE), probability of source resolution, and conditional mean 
squared error (CMSE). Source resolution is defined as the event when both DOAs are estimated 
within one degree of their corresponding true values, i.e., the difference between the true value 
of each DOA and its estimated value is less than 1° x (vr/lSO). The CMSE is defined as 
the expected value of the estimation error conditioned on successful source resolution, i.e.. 



successful source resolution >. The reason for using the CMSE is to fur¬ 


ther investigate the accuracy of the algorithms after making successful detection. We estimate the 
probability of root-swap, subspace leakage, MSE, probability of source resolution, and CMSE 
using the Monte Carlo method with 10® number of trials. Two cases are considered in the 
simulations: 1 ) the two sources are uncorrelated, i.e., r = 0 , and 2 ) the two sources are correlated 
with a correlation coefficient of r = 0.9. 

Eet us start by investigating the probability of root-swap in the root-MUSIC algorithm for 
the case of the uncorrelated sources. The probability of root-swap is estimated using the Monte 
Carlo simulations. Its approximate value is also obtained using (fT5l) . The corresponding curves are 
shown in Eig. [TJ It can be seen that at the low SNR region, the chance that a root-swap occurs is 
quite significant, which results in the performance breakdown of the root-MUSIC algorithm. This 


Febmary 3, 2015 


DRAFT 





16 
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Fig. 1. Probability of root-swap and probability of ML failure versus SNR for uncorrelated sources. 


problem justifies the need for a method to deal with the root-swap phenomenon. In this paper, 
we proposed the root-swap root-MUSIC algorithm which instead of picking the roots closer to 
the unit circle, selects the roots based on the SML criterion. In Fig. [U we thus also draw a 
curve which shows the probability that the selected roots by the ML criterion include a noise 
root. This situation is considered as a breakdown, and therefore, the corresponding probability 
is called the probability of ML failure. As can be seen, this probability is significantly smaller 
than the probability of root-swap. As a result, it is expected that the root-swap root-MUSIC 
algorithm outperforms the conventional root-MUSIC method. This will be shown in the rest of 
this section. 
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The subspace leakage in the two-step root-MUSIC algorithm for the case of the uncorrelated 
sources is investigated next. The expected value of the subspace leakage is estimated using 
(|2^ and the Monte Carlo simulations. The approximate value for the subspace leakage is also 
obtained from the theoretical derivations in (1271) and (l29l) . The value of 7 is fixed at 0.5. The 
results are shown in Fig. [2l The solid lines represent the subspace leakage at the first step, 
and the dashed lines depict the subspace leakage at the second step of the proposed two-step 
root-MUSIC algorithm. It can be seen that the curves obtained from the simulations are very 
close to those obtained from our theoretical derivations at high SNR values. At the low SNR 
region, the curve associated with the theoretical approximation at the second step deviates from 
the curve obtained by simulations. The reason is that in the derivations, the first order Taylor 
series expansion is used. More accurate results can be obtained by using higher order Taylor 
series. However, the computations can become intractable. In Fig. [2l it can be observed from 
both theoretical and simulation results that the subspace leakage from the modified covariance 
matrix at the second step is significantly smaller than the subspace leakage from the sample 
data covariance matrix at the first step. This is achieved by removing the undesirable terms from 
the sample data covariance matrix leading to an estimate of the signal projection matrix that is 
closer to the true signal projection matrix, which is equivalent to a lower subspace leakage at 
the second step. 

We next consider the performance of the proposed two-step algorithm when applied to the root- 
MUSIC dll, unitary root-MUSIC ll22ll . improved unitary root-MUSIC with pseudo-noise resam¬ 
pling lfT2l . root-swap unitary root-MUSIC, and root-swap unitary root-MUSIC with pseudo-noise 
resampling methods. The unitary root-MUSIC algorithm takes benefit from the forward-backward 
averaging [[Tl which is approximately equivalent to doubling the number of samples. For the cases 
that the pseudo-noise resampling is used, P represents the number of times that the resampling 
process has been performed. In the figures, the root-MUSIC, unitary root-MUSIC, and root- 
swap unitary root-MUSIC methods are denoted by R-MUSIC, UR-MUSIC, and RSUR-MUSIC, 
respectively. The value of the scaling factor 7 is obtained by minimizing the SML function as 
described in the two-step root-MUSIC method. In the root-swap algorithm, the parameters p 
and q are set to p = 1 and g = 0, which means the closest root to the unit circle is picked 
up and paired with other roots one at a time in order to find the pair of DOA estimates that 
minimizes the SML function. In this case, the number of different possible combinations of the 
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Fig. 2. Subspace leakage versus SNR for uncorrelated sources. The solid and dashed lines represent the subspace leakage at 
the first and second steps of the proposed two-step root-MUSIC algorithm, respectively. 


roots is Nr = 8. The number of samples used for the pseudo-noise resampling method is set to 
P = 50. Aeeording to our simulations, using more number of samples would not yield in any 
eonsiderable improvement in the performance. 

The MSB versus SNR performance of the methods tested for the case of the uncorrelated 
sources is presented in Fig. [3l The corresponding CRB [|^ is also shown in the figure. For 
the R-MUSIC method, the modification of the covariance matrix in the second step of the 
introduced two-step method shifts the MSB curve by almost half a dB to the left. For the 
UR-MUSIC method the improvement is more significant and is about one dB. For the rest of 
the methods, there is no considerable change in the MSB performance. However, as it will be 
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Fig. 3. MSE versus SNR for uncorrelated sources. The solid and dashed lines are based on the first and second steps of the 
proposed two-step method, respectively. The methods used in the two-step algorithm are R-MUSIC, UR-MUSIC, and RSUR- 
MUSIC methods. P is the number of samples used for the pseudo-noise resampling algorithm. 


shown in the next figures, the modifieation of the eovarianee matrix has benefits in terms of the 
CMSE performanee and probability of souree resolution for these methods. It ean also be seen 
from Fig. [3] that the proposed RSUR-MUSIC algorithm performs about 2 dB better than the UR- 
MUSIC method, while imposing only a small amount of eomputational eomplexity for evaluating 
the SML funetion for Nr = 8 different eombinations of the roots. The best performanee is 
aehieved by the RSUR-MUSIC algorithm eombined with the pseudo-noise resampling method. 

Fig. m shows probability of souree resolution versus SNR for the uneorrelated sourees. For the 
R-MUSIC method, the seeond step of the two-step algorithm improves the performanee by 1 to 
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Fig. 4. Probability of source resolution versus SNR for uncorrelated sources. The solid and dashed lines are based on the first 
and second steps of the proposed two-step method, respectively. The methods used in the two-step algorithm are R-MUSIC, 
UR-MUSIC, and RSUR-MUSIC methods. 


2 dB. The rest of the algorithms have almost the same performanee with the root-swap based 
methods slightly outperforming the other algorithms at low SNR values. It is observed that the 
seeond step of the two-step algorithm results in about 1 dB improvement in the performanee. 

Finally, Fig. [5] illustrates the performanee of the algorithms tested for the uneorrelated sourees 
in terms of the CMSE. The R-MUSIC method is signifieantly improved by the two-step method 
with an improvement ranging from 5 dB at low SNR values to 1 dB at high SNR values. The 
rest of the algorithms show similar performanee, and the applieation of the two-step method 
leads to up to 2 dB improvement in the CMSE performanee. 
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Fig. 5. CMSE versus SNR for uncorrelated sources. The solid and dashed lines are based on the first and second steps 
of the proposed two-step method, respectively. The methods used in the two-step algorithm are R-MUSIC, UR-MUSIC, and 
RSUR-MUSIC methods. 


The results for the ease of the eorrelated sourees with r = 0.9 are depleted in Figs. [6] to 
[TO Similar observations are made from these figures as those discussed for the case of the 
uncorrelated sources. Compared to the uncorrelated case, the performance breakdown occurs at 
a higher SNR value. This makes the importance of the improved methods more significant, as 
there is a higher chance that the actual SNR of a system falls in the breakdown region. As 
seen from the figures for the correlated sources, the proposed methods prove to be helpful in 
dealing with the performance breakdown problem. The gain obtained by the improved methods 
is also more significant compared to the case of the uncorrelated sources. For instance, the MSE 
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Fig. 6. Probability of root-swap and probability of ML failure versus SNR for correlated sources with r = 0.9. 

improvement aehieved by the two-step root-MUSIC method for the uneorrelated sourees is about 
half a dB, while in the ease of the correlated sources, the MSE curve is shifted by more than 
2 dB to the left. Similarly, more significant performance gains are obtained for the probability 
of source resolution and also the CMSE. 

VI. Conclusion 

The performance breakdown of the subspace based DOA estimation methods in the threshold 
region where the SNR and/or sample size is low has been studied in this paper. The subspace 
leakage as the main cause of the performance breakdown was formally defined and theoretically 
derived. The two-step algorithm has been proposed in order to reduce the amount of subspace 
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Fig. 7. Subspace leakage versus SNR for correlated sources with r = 0.9. The solid and dashed lines represent the subspace 
leakage at the first and second steps of the proposed two-step R-MUSIC algorithm, respectively. 


leakage. The introdueed method is based on estimating the DOAs at the first step and modifying 
the covariance matrix using the estimated DOAs at the second step. We have theoretically 
derived the subspace leakage at both steps, and have shown that the subspace leakage is reduced 
at the second step of the proposed method leading to better performance. The algorithm can 
also be extended to the third step by further modifying the covariance matrix based on the 
improved estimates obtained at the second step. We have investigated the performance of the 
algorithm for further steps through simulations (not included in the paper). However, the achieved 
improvement is marginal and does not justify the added complexity. The behavior of the root- 
MUSIC algorithm in the threshold region has been also studied, and a phenomenon called 
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Fig. 8. MSE versus SNR for correlated sources with r = 0.9. The solid and dashed lines are based on the first and second 
steps of the proposed two-step method, respectively. The methods used in the two-step algorithm are R-MUSIC, UR-MUSIC, 
and RSUR-MUSIC methods. 


root-swap has been observed to contribute to the performance breakdown. Then, an improved 
method has been introduced to remedy this problem by considering different combinations of the 
roots and picking up the one that minimizes the SML function. The performance improvement 
achieved by the proposed methods has also been demonstrated using numerical examples and 
simulation results. We also combined the proposed algorithms with the previously introduced 
methods in the literature, which resulted in further improvement in the performance. 
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Fig. 9. Probability of source resolution versus SNR for correlated sources with r = 0.9. The solid and dashed lines are based 
on the first and second steps of the proposed two-step method, respectively. The methods used in the two-step algorithm are 
R-MUSIC, UR-MUSIC, and RSUR-MUSIC methods. 


Appendix A 

Probability of Root-Swap Approximation 

The root-swap is defined as the event when at least one of the estimated signal roots Zk 
0- < k < K) has a smaller magnitude than the magnitude of an estimated noise root Zm 
(K + 1 < m < M — 1), i.e., fk < fm- Let us denote the probability of the event that < Tm by 
Pkm- The eomplement of this event represents the ease when the fc-th estimated signal root has 
not been swapped with the m-th estimated noise root, and its probability is given by 1 —pkm- Let 
us denote the probability of root-swap by P(root-swap). The eomplement of the root-swap event 
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Fig. 10. CMSE versus SNR for the correlated sources with r = 0.9. The solid and dashed lines are based on the first and second 
steps of the proposed two-step method, respectively. The methods used in the two-step algorithm are R-MUSIC), UR-MUSIC, 
and RSUR-MUSIC methods. 


is the event when none of the estimated signal roots has been swapped with an estimated noise 
root, and its probability is given by 1 — P(root-swap). Assuming that the individual root-swap 
events are independent from eaeh other, we have 

K M-l 

1 - P(root-swap) = JJ JJ (1 - pkm) • (31) 

k=l m=K+l 

In the sequel, we derive an approximation for pkm- Noting that = 1 for the true signal 
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roots, we have 


Pkm = P{rm> fk) 

P '^m) • 


(32) 


In order to proeeed with the eomputation of pkm, we eonsider the distributions of Ar^ and 
Arfc. It is shown in ffT9l that Avk (I < k < K) follows the — ((Tfc/\/2) (2(M — K) — 1) 

distribution where (^) denotes a ehi-square distribution with i degrees of freedom and al is 
given by 


err 


CTt, = 


N 


K 


Am- 


M-K+i 


tr - <Tl) 




(33) 


where is the true projection matrix of the noise subspace and is given by (l30l) . 

We next consider the distribution of Ar^- In [|T9l . the distribution of Ar^ is computed using a 
second order Taylor expansion of the estimated root-MUSIC polynomial around the true signal 
roots (which are located on the unit circle). The computation of the distribution of Ar^ requires 
the analysis to be performed around the true noise roots which are located inside the unit circle. 
The second order expansions of a{zk) and a^(i^^) around the true signal root Zk are given by 

m 


a{zk) ^ Ofe + ja^^'^Auk + a^^^Ar^ 

~ af - ja^k^^Auk - a^k^^ Avk (34) 


where Aojk is the difference between the angle of the k-\h estimated root and the angle of the 
corresponding true root. For the m-th noise root, let be defined as 

A [1, • • • , (35) 

where ojm is the angle of Zm- Let also be defined similar to (l30l) with ojk replaced with 
ujjn- Then, the second order expansions of a(Am) and a^{z^) around the true noise root z^ are 
given by 

a{zm) ^ + ja^Aum + ) 

~ - jaP^Aum - (36) 
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where Rm is a M x M diagonal matrix with its diagonal elements equal to 1, •'' > • 

Sinee the Taylor expansion for the steering veetors of the roots on the eirele and the expansion 
for the roots inside the eirele, i.e., (|3^ and (|3^ have similar struetures, it is reasonable to assume 
that Arfc and Ar^/^m also have similar distributions. Then, the varianee of Ar^ is in the order 
of the varianee of Ar^ multiplied by r^. Sinee < 1, the varianee of Ar^ is smaller than the 
varianee of Ar^. In order to simplify the eomputation of pkm, we ignore the effect of Ar^, and 
approximate pkm by 

Pkm ~ P (-Arfc > 1 - r^). (37) 


This is equivalent to using the probability P {fk < Vm) as an approximation for pkm- Since we 
have the distribution of Ar^, we can compute pkm using (l37l) . When M — K 3> 1, Ar^ follows 
approximately a normal distribution AT {^—ak\/M — K — (3/4), lfT9l . Using (l37]) . the 

probability pkm can be approximated by 


Pkn 


Q 


l-Tm- crk\/M - K - (3/4) A 




■ 


(38) 


Finally, the approximation of the probability of root-swap P(root-swap) is found by using the 
approximation (1381) . the expression (l3TI) . and the fact that Q{—x) = 1 — Q{x) as 

■1 + + crk\/ M — K — (3/4) \ 


P (root-swap) 


K M-l 

1-n n e 


k=l m=K-\-l 




V 


(39) 


It completes the derivation. 


Appendix B 

Subspace Leakage at the First Step 

Let us start with the computation of pi. Let AP = P — P he the estimation error of the signal 
projection matrix. Then, using the properties that P^ = P and Tr {P} = K, the expression (1231) 
for the first step of the two-step root-MUSIC algorithm can be written as 

Pi = l-2Tr{(P + AP)P} 

= l-i(A' + Tr{APP}) 

= -7,Tr{APP}. (40) 
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It is shown in IfT^ that the series expansion of P based on AR is given by 

P = P + 6P + --- + S'^P + --- (41) 

where 

6P = P^ARV^ + V^ARP^ (42) 

and the rest of the terms are related by the following reeurrenee 

6^P = -P^ {S^-^P) ARV^ + P^AR (S^-^P) 

-V^AR {S^-^P) P^ + {S^-^P) ARP^ 

n—1 

- ^ P (S^P) {S^-^P) P 

i=l 

n—1 

+ ^ P^ (S^P) (5”-*P) P^. (43) 

i=l 

The following lemma will be further used. 

Lemma 1. The columns ofV^ belong to the signal subspace, i.e., PV^ = Vl 

Proof: The proof follows by multiplying P by and then substituting P with EE^ and 
with dia). ■ 

In a similar way to Lemma [U it ean also be shown that 

yyt (44) 

Using (l40l) . the series expansion of P in (|4TI) . expressions (l42l) and (1431) up to the 5^P term, 

and the faets that PP^ = P^P = 0 and PP = P, we ean write pi as 

p, = -^Tr{-P{6P){6P)}. (45) 

Then, pi is eomputed by substituting (l42l) in (1451) . using P-‘“P-‘- = P-*-, and Lemma [Has 

Pi = ;^Tr|p (P^APV^ + U^^APP^) (P^APV^ + V^ARP^) } 

= ^Tr{PV^APP^P^APVfj 

= ^TrjU^APP^APy^} . (46) 
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Computation of the expected value of the subspace leakage requires considering the statistical 
properties of AR. We use the following two properties in our derivations llT9l . 


Lemma 2. For all matrices Ai, A 2 E we have 

E {ARAiAR} = ^Tr{RAi} R (47) 

and 

E{Tr{ARAi]Tr{ARA 2 ] ] = ^Tr{RAiRA 2 ]. (48) 

Using (l46l) and Wh . the expected value of pi can be computed as 

E{pi} = {V^E [ARP^AR] 

= —Tr I—Tr { RP^ ] RV^ \ 

K { N J 

= j^Tr{P^R}Tr{V^V^R} . (49) 

Since the range space of the matrix A is the same as the signal subspace, we have P^A = 0. 
As a result, Tr can be simplified as 

Tr {P^R} = Tr {ASA^ + a^lM) } 

= Tr{alP^}=alTr{lM-P} 

= al {M-K). (50) 


Furthermore, using (l25l) and the fact that the eigenvectors of R are orthonormal, the product 
V^V^R can be written as 


K 


v^v^r = J2 




M—K+k H 

2^kek 


which results in 


TrjvV^P} = 


{y^M-K+k - crl) 

y^M-K+k 


K 


{y^M-K+k - 0-n) 

Finally, E {pi} is obtained by substituting (l50l) and (15^ in (l49l) as 


E{pi] = 


al (M - K) 


K 


y^M- 


M-K+k 


^ {Xu-K+k - cr^) 


2 ■ 


(51) 


(52) 


(53) 
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Appendix C 

Subspace Leakage at the Second Step 

The subspace leakage at the second step of the two-step root-MUSIC algorithm can be obtained 
through the same steps taken for the computation of pi. Referring to (|4^ . the subspace leakage 
P 2 is given by 

P2 = (54) 

( 2 ) 

where AP^^^ = R — P is the estimation error of the covariance matrix at the second step of 
the algorithm. Using (fT3l) . the estimation error AP^^^ is given by 

AP(2) = AP - 7 (T + T^) . (55) 

Recalling (fTOl) . we have T = PaRPa- 

Consider the first order Taylor series expansion of Pa around the true DOAs given by 


where Pa = A {A^ A) ^ A^ is equal to the true signal projection matrixj, i.e., jt'a 
dP is given by 


Pa^ Pa + dP (56) 

i)ji], i.e., PA = P, and 


K 


dP = 


k=l 


duk 


(57) 


Here Acok ~ djk — ojk is the estimation error of ojk with Afc = 27i{d/\) sm{9k). 

Note that for any square and invertible matrix B, the partial derivative of with respect 
to the variable u is given by 1(241 

= -B-^^B-\ (58) 

OU OOJ 

Using (l58l) . the partial derivative dPA/duk can be computed as 

OPa d{A^A) „ A ( ah A\-^ t'^A 


duk 


dui 


(A^A) A^ + A- 




dujk 

H 


-A^ +A(A^A)’ 


H 


duk 
+A (A^A) 


dA 

duk 


H 


A + A 


H 


duk, 

duk 


A^A)"^A' 


H a\-^ I 

duk 


H 


(59) 


'Note that although Pa is equal to P, the estimates Pa and P are obtained in different ways and are not essentially equal 
to each other. 
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Then, using (l20l) and P = A A) ^ A^, the partial derivative dPA/d^k is given by 


dP, 

dojk 


= P 


dA 

duk 


{A^A) ^A^ + A{A^A) 


H A\-^ 


dA 

duk 


H 


The estimation error of Uk, i.e., Auk in (1571) . ean be written based on AR as [[T^ 

a[^^^P^ARV^ak - ARP^a^f^^ 


Auk = 


2 ] I 


The first order Taylor series expansion of is obtained using (fT2l) and (15^ as 


(60) 


(61) 


P^^Pi-dP 


(62) 


where P\ = Im - Pa- 

The matrix T can be then computed using expressions (fTOl) . (15^ . and (l62l) with keeping only 
the first order terms and noting that Pa = P, Pa = P^^ PRP^ = 0 as 


T = {Pa + dP) {R + AR) {P\ - dP) 
^ -PRdP + PARP^ + dPRP^. 


(63) 


We can now compute p 2 using expressions (l54l) . (1551) . and (1631) as 

P2 = ;^Tr|Fi (AP- 7 (T + T^))P^ (AP - 7 (T + T^)) 

= (AP - 7 ( - PRdP + PARP^ + dPRP^ - dPRP + P^ARP 

+P^RdP))P^{AR - 7( - PRdP + PARP^ + dPRP^ 

-dPRP + P^ARP + P^RdP))V^y (64) 

Then, using expressions (I57l) . (l60l) . and the fact that PP^ = P^P = P'^ = P'^V^ = 0 to 

eliminate the terms that equal zero, p 2 is computed as 

P2 = -l{- PRdP + PARP^ + dPRP^)) 

X P^ (AP - 7 ( - dPRP + P^APP + P^PdP)) I. (65) 
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Expanding the terms in (1651) and using the faet that PV^ = P = results in the following 
expression for p 2 

P2 = ^1:x{v^ARP^ARV^ --fi-V^ARP^dPRV^ + V^ARP^^RV^ 
+V^/\RP^RdPV^ - V^RdPP^/\RV^ + V^^RP^^RV^ 
+V^dPRP^^RV^) +-i‘^{V^RdPP^dPRV^ - V^RdPP^^RV^ 
-V^RdPP^RdPV^ - V^^RP^dPRV^ + V^ARP^^RV^ 
+V^^RP^RdPV^ - VUPRP^dPRV^ + VUPRP^^RV^ 
+VUPRP^RdPV^)y 

(66) 

By reordering the terms in (1^ . the subspaee leakage p 2 can be further rewritten as 

P 2 = ^Tr| (1 - 27 + 7 ^) V^ARP^ARV^ + ( 7 ^ - 7 ) ( - V^ARP^dPRV^ 
+V^ARP^RdPV^ - V^RdPP^ARV^ + vUpRP^ARV^) 

RdPP^dPRV^ - V^RdPP^RdPV^ - VUPRP^dPRV^ 
+VUPRP^RdPV^)y (67) 

The terms multiplied by ( 7 ^ — 7 ) in (1671) can be simplified using expressions (l24l) . (@41), and 
the fact that P^V = 0 as 

-V^ARP^dP {V + allM) + V^ARP^ {V + allu) dPV^ 

-V^ {V + (tIIm) dPP^ARV^ + V^dP [V + ct^Jm) P^ARV^ 

= -V^ARP^dPP - PdPP^ARV^ (68) 

In a similar way, the terms multiplied by 7 ^ in (l67l) can be simplified to 

V^RdPP^dP {V + allu) - V^RdPP^ {V + ct^Jm) dPV^ 
-V^dPRP^dP {V + allu) + VUPRP^ [V + (jIIm) dPV^ 

= V^RdPP^dPP - V^dPRP^dPP 
= V^ {V + allu) dPP^dPP - VUP {V + allM) P^dPP 
= PdPP^dPP (69) 
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which using the fact that P^dPP^ = 0 (see (1571) and (l60l) ) can be further simplified to 
PdPP^dPP = {Im - P^) dPP^dP {Im - P^) 

= dPP^dP. ( 70 ) 

Finally, using expressions (l46l) . (1671) . (l68l) . (170l) . and Lemma HI the subspace leakage p 2 is 
computed as 

P2= ( 1 -27 + 72)^1+ ^il^^i?e{Tr{ VtAi^P^dP}} + ^Tr{dPP^dP}. ( 71 ) 

Computation of the expected value of p 2 involves finding the expected value of the two trace 
functions in (ITTl) . Using expressions (1571) and (l60l) . the expected value of the first trace function 
in (ITTl) is given by 

K 


dA 


E {Tr{V^ARP^dP}} = P <j Tr <1 AP^ P^^ (A^A) AukV 




k=l 


(72) 


Then, by substituting (l6Tl) in (l72l) . we have 


K 


dA 


E {Tr{V^ARP^dP}} = P<j Tr<' ^APP^^ (A^A) A^V^ 


k=l 


X- 


2j ( 


a^dHp±AjiY\f^^ _ a^V^ARP^a^,^^ 


(73) 


The order of the summation and trace operator in (1731) can be swaped. Moreover, the last two 
terms can be written using the trace operator as 

1 


K 


E {Tr{V^ARP^dP}} = eI ^ — 


“* n- I jd-L W 

k=i 2j (rifc P 


X 


TrjAPP^^ (A^A) ^ A^V^I (^Tr |aPUW4^^^^^} 


-Tr \ ARP^a!'^'^ 


(74) 
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The expression in (1741) ean be computed using (l48l) as 


k=l [O'k ^ 


X I Tr<{ {A^A) ^ A^V^RV^akal^'^^P^ 


-Tr<^ RP 


dA 




— (A"A)- A«V^RP 


The second trace function in (fTSl) equals zero as V^RP^ = 0. Then, expression 
rewritten as 


E {Tr{V^ARP^dP}} = 


K ^A^V^RV^ak 


k=l 
r2 D-L 




(75) 
can be 

(76) 


where we used the equality P R = cr^P ■ 

In a similar way, using expressions (1571) and (l60l) . the expected value of the second trace 
function in (ITT]) is given by 


E {Tr{dPP^dP}} = 


K K 


y;y;A(A"A)-' 

( ( k=l i=l 

Then, by substituting (IMl) in (1771) . we have 

K K 


dA 

duk 


H 


(A"A)-'A"Aa.tAa.i 

dui ^ 


(77) 


E{Tr{dPP^dP}} = E<j Tr<j ^^A(A^A) ^ 

k=l i=l 
1 


dA 

duk 


H 


IA« A)-'A" 

doji '' 


X- 


2j (a«"P^o« 


X 


2j (aW"P^af > 


X 


X 


Tr I APV^afc^^^^P^} - Tr | APP^a^^^a^F^ | 
Tr I APV^a.af - Tr | APP^af ^af | 


(78) 
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which is computed using ( 1 ^ and the fact that P^RV^ = 0 as 


K K 


E{Tv{dPPUP}}^^J2J2 






X 


Re |af . (79) 


Finally, the expected value of p 2 for a fixed value of 7 is obtained using expressions (TtT]) . 
(17^ . and (1791) as 


E { p2 } = (1 - 27 + 7^) E {pi} 


2 (7 - 7 ) 
NK 


^ (A^A) ^A^V^RV^au 




k=l 


2j ( aP^P^a^P 


(S^'] 


H 


(A^A) 


-1 




Re |af . 

(80) 


It concludes the derivation. 
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